In [64]:
%matplotlib inline
import matplotlib.pyplot as plt
from datetime import timedelta
import pandas

In [101]:
import sunpy.lightcurve
from sunpy.time import TimeRange, parse_time

In [113]:
launch_date = parse_time('2015/09/15 00:00')

In [12]:
tr = TimeRange('2001/01/01', '2020/01/01')
noaa = sunpy.lightcurve.NOAAIndicesLightCurve.create(tr)
noaa_predict = sunpy.lightcurve.NOAAPredictIndicesLightCurve.create(tr)

In [154]:
l = noaa_predict.data['sunspot']
nl = l.resample('D')
nl = nl.interpolate(method='linear')
nl.plot()
snum = nl[launch_date]
l[l==snum].index.tolist()


Out[154]:
[]

In [40]:
from sunpy.net import hek
client = hek.HEKClient()
result = client.query(hek.attrs.Time(tstart, tend), hek.attrs.EventType('FL'), hek.attrs.FRM.Name=='SSW Latest Events')

In [41]:
goes_class = [event.get('fl_goescls')[0] for event in result]

In [42]:
goes_class.count('M')


Out[42]:
2

In [75]:
def numflares_per_interval(time_range, dt):
    trs = time_range.window(dt, dt)
    classes = ['C', 'M', 'X']
    CClass = []
    MClass = []
    XClass = []
    time = []
    for tr in trs:
        hek_result = client.query(hek.attrs.Time(tr.start, tr.end), hek.attrs.EventType('FL'), hek.attrs.FRM.Name=='SSW Latest Events')
        goes_class_list = [event.get('fl_goescls')[0] for event in hek_result]
        CClass.append(goes_class_list.count('C'))
        MClass.append(goes_class_list.count('M'))
        XClass.append(goes_class_list.count('X'))
        time.append(tr.start)
    data = {"X": XClass, "M": MClass, "C": CClass}
    return pandas.DataFrame(data, index=time)

In [80]:
result = numflares_per_interval(TimeRange('2011/01/01 00:00:00', '2011/06/20 00:00:00'), timedelta(1))

In [82]:
result.plot()


Out[82]:
<matplotlib.axes.AxesSubplot at 0x10d2b97d0>

In [159]:
def cross(series, cross=0, direction='cross'):
    """
    Given a Series returns all the index values where the data values equal 
    the 'cross' value. 

    Direction can be 'rising' (for rising edge), 'falling' (for only falling 
    edge), or 'cross' for both edges
    from http://stackoverflow.com/questions/10475488/calculating-crossing-intercept-points-of-a-series-or-dataframe
    """
    # Find if values are above or bellow yvalue crossing:
    above=series.values > cross
    below=np.logical_not(above)
    left_shifted_above = above[1:]
    left_shifted_below = below[1:]
    x_crossings = []
    # Find indexes on left side of crossing point
    if direction == 'rising':
        idxs = (left_shifted_above & below[0:-1]).nonzero()[0]
    elif direction == 'falling':
        idxs = (left_shifted_below & above[0:-1]).nonzero()[0]
    else:
        rising = left_shifted_above & below[0:-1]
        falling = left_shifted_below & above[0:-1]
        idxs = (rising | falling).nonzero()[0]

    # Calculate x crossings with interpolation using formula for a line:
    x1 = series.index.values[idxs]
    x2 = series.index.values[idxs+1]
    y1 = series.values[idxs]
    y2 = series.values[idxs+1]
    x_crossings = (cross-y1)*(x2-x1)/(y2-y1) + x1

    return x_crossings

In [170]:
solutions = cross(noaa.data['sunspot SWO smooth'],cross=snum)
pandas.Timestamp(solutions[0])


Out[170]:
Timestamp('2004-12-19 11:05:35.999999995')

In [225]:
dt = timedelta(30)
fig = plt.figure(figsize=(20, 15))
fig.add_subplot(1, 1, 1)
noaa.data['sunspot SWO'].plot()
noaa.data['sunspot SWO smooth'].plot()
noaa_predict.data['sunspot'].plot()
plt.axvline(launch_date, color='black', linestyle='--', label='Launch date')
plt.legend()
plt.axhline(snum)
for sol in solutions[0:-1]:
    plt.axvline(pandas.Timestamp(sol))
    plt.axvline(pandas.Timestamp(sol) - dt, linestyle='--')
    plt.axvline(pandas.Timestamp(sol) + dt, linestyle='--')
plt.show()



In [226]:
data = numflares_per_interval(TimeRange(pandas.Timestamp(solutions[1])-dt, pandas.Timestamp(solutions[1])+dt), timedelta(1))
data.hist()


Out[226]:
array([[<matplotlib.axes.AxesSubplot object at 0x10dbaf0d0>,
        <matplotlib.axes.AxesSubplot object at 0x1103064d0>],
       [<matplotlib.axes.AxesSubplot object at 0x10db28690>,
        <matplotlib.axes.AxesSubplot object at 0x110319c10>]], dtype=object)

In [227]:
for goes_class in data.columns:
    lc = data[goes_class]
    percent = len(lc[lc == 0]) / np.float(len(lc)) * 100
    print('Percent of days with no ' + goes_class + ' ' + str(percent))


Percent of days with no C 18.3333333333
Percent of days with no M 71.6666666667
Percent of days with no X 96.6666666667

No data before 2011!!!


In [200]:
print(numflares_per_interval(TimeRange('2004/01/01', '2004/02/01'), timedelta(1)))


            C  M  X
2004-01-01  0  0  0
2004-01-02  0  0  0
2004-01-03  0  0  0
2004-01-04  0  0  0
2004-01-05  0  0  0
2004-01-06  0  0  0
2004-01-07  0  0  0
2004-01-08  0  0  0
2004-01-09  0  0  0
2004-01-10  0  0  0
2004-01-11  0  0  0
2004-01-12  0  0  0
2004-01-13  0  0  0
2004-01-14  0  0  0
2004-01-15  0  0  0
2004-01-16  0  0  0
2004-01-17  0  0  0
2004-01-18  0  0  0
2004-01-19  0  0  0
2004-01-20  0  0  0
2004-01-21  0  0  0
2004-01-22  0  0  0
2004-01-23  0  0  0
2004-01-24  0  0  0
2004-01-25  0  0  0
2004-01-26  0  0  0
2004-01-27  0  0  0
2004-01-28  0  0  0
2004-01-29  0  0  0
2004-01-30  0  0  0
2004-01-31  0  0  0

In [ ]: